/*
 * Copyright 2003 Program of Computer Graphics, Cornell University
 *     580 Rhodes Hall
 *     Cornell University
 *     Ithaca NY 14853
 * Web: http://www.graphics.cornell.edu/
 * 
 * Not for commercial use. Do not redistribute without permission.
 */
#include "stdafx.h"

#include "KDTree.h"
#include "PrimitiveList.h"
#include "Ray.h"
#include "Scene.h"
#include "Triangle.h"
#include <float.h> //    FLT_MIN, FLT_MAX, FLT_EPSILON
#include <iostream>

#define KDNGC_HISTORY 4
class KDNodeGroupCombiner {
  int nextReplace;
  KDNodeGroup *putBack;
  struct { KDNodeGroup *group; int size; } history[KDNGC_HISTORY];

  void combineGroups(KDNodeGroup *dest, int destSize, 
    KDNodeHandle **srcRootHandle, KDNodeGroup *src, int srcSize) {
    assert(destSize+srcSize<=static_cast<int>(KDNODEGROUP_COUNT));
    //copy nodes from src to extra space in dest, updating pointers as necessary
    for(int i=0;i<srcSize;i++) {
      KDNode *srcNode = &src->node[i];
      KDNode *destNode = &dest->node[i+destSize];
      destNode->splitValue = srcNode->splitValue;
      int leftType = getNodeType(srcNode->left);
      if (leftType == KDST_LEAF) {
        destNode->left = srcNode->left;
      } else {
        LONG_PTR offset = getNodePtr(srcNode->left)-srcNode;
        assert((i+offset<srcSize)&&(offset>=0));
        destNode->left = makeKDNodeHandle(destNode+offset,leftType);
      }
      int rightType = getNodeType(srcNode->right);
      if (rightType == KDST_LEAF) {
        destNode->right = srcNode->right;
      } else {
        LONG_PTR offset = getNodePtr(srcNode->right)-srcNode;
        assert((i+offset<srcSize)&&(offset>=0));
        destNode->right = makeKDNodeHandle(destNode+offset,rightType);
      }
    }
    //reset sources root to point to the subtree's new location in dest
    int rootType = getNodeType(*srcRootHandle);
    *srcRootHandle = makeKDNodeHandle(&dest->node[destSize],rootType);
  }
public:
  KDNodeGroupCombiner() {
    putBack = NULL;
    for(int i=0;i<KDNGC_HISTORY;i++) {
      history[i].group = NULL; history[i].size = KDNODEGROUP_COUNT;
    }
    nextReplace = 0;
  }

  KDNodeGroup *getNewGroup(PoolAllocator *pool) {
    KDNodeGroup *retval;
    if (putBack) { 
      retval = putBack;
      putBack = NULL;
    } else {
      retval = (KDNodeGroup*)pool->allocate(sizeof(KDNodeGroup),KDNODEGROUP_ALIGNMENT);
    }
    return retval;
  }
  void handleUnderfullGroup(KDNodeHandle **rootHandle, KDNodeGroup *group, int size) {
    //try to find another group to combine this group with 
    int bestIndex = -1;
    int i = nextReplace;
    do {   //find most recent fit
      i = (i-1)&(KDNGC_HISTORY-1);
      if (history[i].size + size < static_cast<int>(KDNODEGROUP_COUNT)) {
        bestIndex = i;
        break;
      }
    } while (i!=nextReplace);
    /*
    int bestIndex = -1;
    int bestSpace = KDNODEGROUP_COUNT;
    for(int i=0;i<KDNGC_HISTORY;i++) { //best fit
      int spareSpace = KDNODEGROUP_COUNT - size - history[i].size;
      if ((spareSpace>=0)&&(spareSpace<bestSpace)) {
        bestIndex = i;
        bestSpace = spareSpace;
      }
    }
    */
    if (bestIndex>=0) {
      combineGroups(history[bestIndex].group,history[bestIndex].size,
        rootHandle,group,size);
      history[bestIndex].size += size;
      assert(putBack==NULL);
      putBack = group;
    } else {
      //otherwise put this group in the history list using round robin replacement
      history[nextReplace].group = group;
      history[nextReplace].size = size;
      nextReplace = (nextReplace<KDNGC_HISTORY-1)?(nextReplace+1):0;
    }
  }
};


KDTree::KDTree(KDBuilder *buildStrategy, Scene *scene)
:nodePool(POOL_SEGMENT_SIZE,sizeof(float))
,primitiveListPool(POOL_SEGMENT_SIZE,SIMD_SIZE)
,fatNodePool(POOL_SEGMENT_SIZE,sizeof(float))
{
  //establish the bounding box
  bounds.resetToEmpty();
  GeometrySet *geom = scene->getAllGeometrySet();
  geom->expandBoxToInclude(&bounds);
  epsBounds = bounds;
//  epsBounds.epsilonExpand(EPSM_BUILD);
  epsBounds.epsilonExpand(EPSM_BIG);
//  buildNode(&rootPtr,buildStrategy,0,bounds,geom);
  buildNode(&preliminaryRootPtr,buildStrategy,0,bounds,*geom);
  //uncompiled geometry is no longer needed at this point
  scene->deleteUncompiledGeometry();
  combiner = new KDNodeGroupCombiner();
//  compileNodesSimplistic(&rootHandle,preliminaryRootPtr);
//  compileNodesChunkedBreadthFirst(&rootHandle,preliminaryRootPtr);
  compileNodesChunkedMostLikely(&rootHandle,preliminaryRootPtr);
  delete combiner;
  combiner = NULL;
  fatNodePool.deallocateAll();
  preliminaryRootPtr = NULL;
  dummy.splitValue = 0.0f;
  dummy.left = dummy.right = NULL;
}

KDTree::~KDTree() 
{
  if (combiner) delete combiner;
}


#define DEBUGCODE(a) a
//#define DEBUGCODE(a)  
void KDTree::buildNode(KDFatNode **res, KDBuilder *builder, int depth, const AxisAlignedBox &bounds, GeometrySet &geom)
{
  KDSplitCode code;
  if (depth >= KDTREE_MAX_DEPTH) {
    printf("Warning tree max depth reached during build\n");
    code = KDSplitCode(KDST_LEAF);
  } else {
    code = builder->split(depth,bounds,geom);
  }
  int type = code.getSplitType();
  if (type==KDST_LEAF) {
    *res = (KDFatNode*)geom.compileToPrimitiveList(&primitiveListPool);
  } else {
    DEBUGCODE(int predLeftTri = builder->buildPredLeftTri);
    DEBUGCODE(int predRightTri = builder->buildPredRightTri);
    DEBUGCODE(int predLeftOther = builder->buildPredLeftOther);
    DEBUGCODE(int predRightOther = builder->buildPredRightOther);
    float splitValue = code.getSplitValue();
    KDFatNode *node = (KDFatNode *)fatNodePool.allocate(sizeof(KDFatNode));
    node->set(code,depth,bounds);
    *res = node;
    AxisAlignedBox childBox = bounds;
    childBox.max[type] = splitValue;
    GeometrySet *childSet = new GeometrySet(geom,false);
    childSet->setSubsetIntersectsBox(geom,childBox);
    DEBUGCODE(int leftTri = childSet->numTriangles());
    DEBUGCODE(int leftOther = childSet->numOthers());
    buildNode(&node->left,builder,depth+1,childBox,*childSet);
    childBox = bounds;
    childBox.min[type] = splitValue;
    childSet->setSubsetIntersectsBox(geom,childBox);
    DEBUGCODE(int rightTri = childSet->numTriangles());
    DEBUGCODE(int rightOther = childSet->numOthers());
    DEBUGCODE(if ((predLeftTri>=0)&&((predLeftTri<leftTri)||(predRightTri<rightTri))) )
    DEBUGCODE(  printf("tri split %d %d did not match prediction %d %d\n",leftTri,rightTri,predLeftTri,predRightTri); )
    DEBUGCODE(if ((predLeftOther>=0)&&((predLeftOther<leftOther)||(predRightOther<rightOther))) )
    DEBUGCODE(  printf("other split %d %d did not match prediction %d %d\n",leftOther,rightOther,predLeftOther,predRightOther); )
    buildNode(&node->right,builder,depth+1,childBox,*childSet);
    delete childSet;
  }
}
#undef DEBUGCODE

void KDTree::compileNodesSimplistic(KDNodeHandle **handle, KDFatNode *fat) 
{
  if (fat==NULL) {
    *handle = makeKDNodeLeafHandle(NULL);
  } else {
    int type = fat->getSplitType();
    if (type == KDST_LEAF) {
      //in this case fat is actually a PrimitiveList pointer
      *handle = makeKDNodeLeafHandle((PrimitiveList *)fat);
    } else {
      KDNode *node = (KDNode *)nodePool.allocate(sizeof(KDNode));
      node->splitValue = fat->splitValue;
      *handle = makeKDNodeHandle(node,type);
      compileNodesSimplistic(&node->left,fat->left);
      compileNodesSimplistic(&node->right,fat->right);
    }
  }
}

void KDTree::compileNodesChunkedBreadthFirst(KDNodeHandle **rootHandle, KDFatNode *rootFat)
{
  //if it is a leaf, just embed it directly, no need to allocate any nodes
  if ((rootFat==NULL)||(rootFat->getSplitType() == KDST_LEAF)) {
    *rootHandle = makeKDNodeLeafHandle((PrimitiveList*)rootFat); 
    return;
  }
  const int arraySize = KDNODEGROUP_COUNT+2;
  //data is a circular array used to hold nodes to be processed in breadth first order
  struct { KDNodeHandle **handle; KDFatNode *fat;} data[arraySize];
  int cur = 0;
  int limit = 1;
  data[cur].handle = rootHandle;
  data[cur].fat = rootFat;
  KDNodeGroup *grp = combiner->getNewGroup(&nodePool);
  int i=0;
  for(;i<static_cast<int>(KDNODEGROUP_COUNT);i++) {
    if (cur==limit) break; //if no more nodes to process then break out of loop
    KDNodeHandle **handle = data[cur].handle;
    KDFatNode *fat = data[cur].fat;
    assert(fat->getSplitType()!=KDST_LEAF);
    KDNode *node = &grp->node[i];
    node->splitValue = fat->splitValue;
    *handle = makeKDNodeHandle(node,fat->getSplitType());
    cur = (cur<arraySize-1)?(cur+1):0;
    //process left child, if its a leaf just embed it directly
    if ((fat->left==NULL)||(fat->left->getSplitType()==KDST_LEAF)) {
      node->left = makeKDNodeLeafHandle((PrimitiveList*)fat->left);
    } else { //otherwise enqueue it for later processing in breadth first order
      data[limit].handle = &node->left;
      data[limit].fat = fat->left;
      limit = (limit<arraySize-1)?(limit+1):0;
      assert(limit!=cur);
    }
    //process right child, if its a leaf just embed it directly
    if ((fat->right==NULL)||(fat->right->getSplitType()==KDST_LEAF)) {
      node->right = makeKDNodeLeafHandle((PrimitiveList*)fat->right);
    } else { //otherwise enqueue it for later processing in breadth first order
      data[limit].handle = &node->right;
      data[limit].fat = fat->right;
      limit = (limit<arraySize-1)?(limit+1):0;
      assert(limit!=cur);
    }
  }
  if (i<static_cast<int>(KDNODEGROUP_COUNT)) { //try to combine underfull groups to save memory
    combiner->handleUnderfullGroup(rootHandle,grp,i);
  } else {
    //after packing as many nodes as will fit into group, the recurse to process any remaining
    //nodes
    while(cur!=limit) {
      compileNodesChunkedBreadthFirst(data[cur].handle,data[cur].fat);
      cur = (cur<arraySize-1)?(cur+1):0;
    }
  }
}

void KDTree::compileNodesChunkedMostLikely(KDNodeHandle **rootHandle, KDFatNode *rootFat)
{
  //if it is a leaf, just embed it directly, no need to allocate any nodes
  if ((rootFat==NULL)||(rootFat->getSplitType() == KDST_LEAF)) {
    *rootHandle = makeKDNodeLeafHandle((PrimitiveList*)rootFat); 
    return;
  }
  const int arraySize = KDNODEGROUP_COUNT+1;
  //data is a circular array used to hold nodes to be processed in breadth first order
  struct { KDNodeHandle **handle; KDFatNode *fat;} data[arraySize];
  int num = 1;
  data[0].handle = rootHandle;
  data[0].fat = rootFat;
  KDNodeGroup *grp = combiner->getNewGroup(&nodePool);
  int i=0;
  for(;i<static_cast<int>(KDNODEGROUP_COUNT);i++) {
    if (num==0) break; //if no more nodes to process then break out of loop
    //search for node with largest surface area
    F32vec1 biggestArea(data[0].fat->area);
    int biggestIndex = 0;
    for(int j=1;j<num;j++) {
      F32vec1 testArea(data[j].fat->area);
      if (_mm_ucomigt_ss(testArea,biggestArea))
        biggestIndex = j;
      biggestArea = simd_max(biggestArea,testArea);
    }
    //process it, then remove it by swapping it with last element in the area
    KDNodeHandle **handle = data[biggestIndex].handle;
    KDFatNode *fat = data[biggestIndex].fat;
    assert(fat->getSplitType()!=KDST_LEAF);
    num--;
    data[biggestIndex] = data[num];
    KDNode *node = &grp->node[i];
    node->splitValue = fat->splitValue;
    *handle = makeKDNodeHandle(node,fat->getSplitType());
    //process left child, if its a leaf just embed it directly
    if ((fat->left==NULL)||(fat->left->getSplitType()==KDST_LEAF)) {
      node->left = makeKDNodeLeafHandle((PrimitiveList*)fat->left);
    } else { //otherwise enqueue it for later processing in breadth first order
      data[num].handle = &node->left;
      data[num].fat = fat->left;
      num++;
      assert(num<=arraySize);
    }
    //process right child, if its a leaf just embed it directly
    if ((fat->right==NULL)||(fat->right->getSplitType()==KDST_LEAF)) {
      node->right = makeKDNodeLeafHandle((PrimitiveList*)fat->right);
    } else { //otherwise enqueue it for later processing in breadth first order
      data[num].handle = &node->right;
      data[num].fat = fat->right;
      num++;
      assert(num<=arraySize);
    }
  }
  if (i<static_cast<int>(KDNODEGROUP_COUNT)) { //try to combine underfull groups to save memory
    combiner->handleUnderfullGroup(rootHandle,grp,i);
  } else {
    //after packing as many nodes as will fit into group, the recurse to process any remaining
    //nodes
    for(int k=0;k<num;k++) {
      compileNodesChunkedMostLikely(data[k].handle,data[k].fat);
    }
  }
}


class KDTStats {
public:
  int numNodes, numLeafNodes, numEmptyLeafNodes;
  int sumLeafDepth, minLeafDepth, maxLeafDepth;
  int sumLeafSize, maxLeafSize;
  KDTStats() {
    numNodes = numLeafNodes = numEmptyLeafNodes = 0;
    sumLeafDepth = maxLeafDepth = 0;
    minLeafDepth = KDTREE_MAX_DEPTH;
    sumLeafSize = maxLeafSize = 0;
  }
  void collectStats(KDNodeHandle *handle, int depth);
};

void KDTStats::collectStats(KDNodeHandle *handle, int depth)
{
  int type = getNodeType(handle);
  KDNode *node = getNodePtr(handle);
  numNodes += 1;
  if (type == KDST_LEAF) {
    sumLeafDepth += depth;
    minLeafDepth = minv(minLeafDepth,depth);
    maxLeafDepth = maxv(maxLeafDepth,depth);
    int size = 0;
    for(PrimitiveList *list=(PrimitiveList*)node; list!=NULL; list=list->next) {
      if (list->isQuadPackedTriangleList()) {
        size += ((QuadPackedTriangleList*)list)->computeLength();
        break;  //always last in list and may not have a next pointer if USE_COMPILED_SHARED_TRIANGLE is set
      } else {
        size += list->length;
      }
    }
    sumLeafSize += size;
    maxLeafSize = maxv(maxLeafSize,size);
    if (size==0) numEmptyLeafNodes += 1;
    numLeafNodes += 1;
  } else {
    collectStats(node->left,depth+1);
    collectStats(node->right,depth+1);
  }
}

void KDTree::printTreeStatistics() 
{
  KDTStats stats;
  stats.collectStats(rootHandle, 0);
  printf("tree contains %d nodes, nodes %.3g MB, primitives %.3g MB\n",stats.numNodes,
    nodePool.getTotalSize()/(1024.0*1024.0),
    primitiveListPool.getTotalSize()/(1024.0*1024.0));
  printf("  %d leaf nodes (%.3g%% empty), depth min %d, max %d, avg %.3g\n",
    stats.numLeafNodes,stats.numEmptyLeafNodes*100.0/stats.numLeafNodes,
    stats.minLeafDepth, stats.maxLeafDepth, stats.sumLeafDepth/(double)stats.numLeafNodes);
  printf("  leaf size avg %.3g, non-empty avg %.3g, max %d\n",
    stats.sumLeafSize/(double)stats.numLeafNodes,
    stats.sumLeafSize/(double)(stats.numLeafNodes-stats.numEmptyLeafNodes),
    stats.maxLeafSize);
  fflush(stdout);
}

KDSplitCode KDBFlatList::split(int depth, const AxisAlignedBox &bounds, const GeometrySet &geom)
{
  return KDSplitCode(KDST_LEAF); //we only create leaf nodes for a flat list
}

KDSplitCode KDBNaiveMiddleSplit::split(int depth, const AxisAlignedBox &bounds, const GeometrySet &geom)
{
  if ((depth>=maxDepth)||(geom.size()==0)) {
    return KDSplitCode(KDST_LEAF); //we only create leaf nodes for a flat list
  } else { //split axis with the largest extent
    float dx = bounds.max[0] - bounds.min[0];
    float dy = bounds.max[1] - bounds.min[1];
    float dz = bounds.max[2] - bounds.min[2];
    int type;
    if ((dx > dy)&&(dx > dz)) type = KDST_SPLIT_X;
    else if (dy > dz) type = KDST_SPLIT_Y;
    else type = KDST_SPLIT_Z;
    return KDSplitCode(type,0.5f*(bounds.min[type]+bounds.max[type]));
  }
}

KDSplitCode KDBNaiveMedianSplit::split(int depth, const AxisAlignedBox &bounds, const GeometrySet &geom)
{
  if (geom.size()<=maxElements) {
    return KDST_LEAF; //we only create leaf nodes for a flat list
  } else { //split axis with the largest extent
    float dx = bounds.max[0] - bounds.min[0];
    float dy = bounds.max[1] - bounds.min[1];
    float dz = bounds.max[2] - bounds.min[2];
    int type;
    if ((dx > dy)&&(dx > dz)) type = KDST_SPLIT_X;
    else if (dy > dz) type = KDST_SPLIT_Y;
    else type = KDST_SPLIT_Z;
    
    KDFixedSplitEstimator splits(bounds);
    geom.estimateFixedSplits(splits);
    splits.integrateCounts();
    int totalSize = splits.getTotalTri() + splits.getTotalOther();
    int bestIndex = 32;
    int bestVal,i;
    switch(type) {
    case KDST_SPLIT_X:
      bestVal = maxv(splits.getLeftTriCountX(32)+splits.getLeftOtherCountX(32),
                        splits.getRightTriCountX(32)+splits.getRightOtherCountX(32));
      for(i=1;i<32;i++) {
        int testVal = maxv(splits.getLeftTriCountX(32-i)+splits.getLeftOtherCountX(32-i),
                          splits.getRightTriCountX(32-i)+splits.getRightOtherCountX(32-i));
        if (testVal < bestVal) {
          bestIndex = 32-i; bestVal = testVal;
        }
        testVal = maxv(splits.getLeftTriCountX(32+i)+splits.getLeftOtherCountX(32+i),
                      splits.getRightTriCountX(32+i)+splits.getRightOtherCountX(32+i));
        if (testVal < bestVal) {
          bestIndex = 32+i; bestVal = testVal;
        }
      }
      break;
    case KDST_SPLIT_Y:
      bestVal = maxv(splits.getLeftTriCountY(32)+splits.getLeftOtherCountY(32),
                        splits.getRightTriCountY(32)+splits.getRightOtherCountY(32));
      for(i=1;i<32;i++) {
        int testVal = maxv(splits.getLeftTriCountY(32-i)+splits.getLeftOtherCountY(32-i),
                          splits.getRightTriCountY(32-i)+splits.getRightOtherCountY(32-i));
        if (testVal < bestVal) {
          bestIndex = 32-i; bestVal = testVal;
        }
        testVal = maxv(splits.getLeftTriCountY(32+i)+splits.getLeftOtherCountY(32+i),
                      splits.getRightTriCountY(32+i)+splits.getRightOtherCountY(32+i));
        if (testVal < bestVal) {
          bestIndex = 32+i; bestVal = testVal;
        }
      }
      break;
    case KDST_SPLIT_Z:
      bestVal = maxv(splits.getLeftTriCountZ(32)+splits.getLeftOtherCountZ(32),
                      splits.getRightTriCountZ(32)+splits.getRightOtherCountZ(32));
      for(i=1;i<32;i++) {
        int testVal = maxv(splits.getLeftTriCountZ(32-i)+splits.getLeftOtherCountZ(32-i),
                          splits.getRightTriCountZ(32-i)+splits.getRightOtherCountZ(32-i));
        if (testVal < bestVal) {
          bestIndex = 32-i; bestVal = testVal;
        }
        testVal = maxv(splits.getLeftTriCountZ(32+i)+splits.getLeftOtherCountZ(32+i),
                      splits.getRightTriCountZ(32+i)+splits.getRightOtherCountZ(32+i));
        if (testVal < bestVal) {
          bestIndex = 32+i; bestVal = testVal;
        }
      }
      break;
    }
    if (bestVal==totalSize) {//didn't manage to reduce list
      return KDSplitCode(KDST_LEAF);
    }
    return splits.getSplitCode(type,bestIndex);
  }
}

KDBGreedyCostSplit::KDBGreedyCostSplit() {
  for(int i=0;i<KDFIXEDSPLIT_NUM;i++) {
    relSplitTable[i] = float(i)/KDFIXEDSPLIT_NUM;
  }
}
static inline float nodeCost(int numTri, int numOther) {
  int numTriPack = (numTri+3)&~0x03;
  return static_cast<float>( numTriPack+numOther );
}
static inline float splitCost(float leftProb, int leftTri, int leftOther,
                              float rightProb, int rightTri, int rightOther,
                              int totalTri, int totalOther) {
//  const float splitNodeCost = 2.0f;  //base cost for splitting node
//  const float duplicateCost = 0.25f; //extra cost for duplicating primitives in both children
//  const float reductionBenefit = 0.05f; //extra benefit for reducing max primitive count
  const float splitNodeCost = 0.75f;  //base cost for splitting node
  const float duplicateCost = 0.15f; //extra cost for duplicating primitives in both children
//  const float reductionBenefit = 0.05f; //extra benefit for reducing max primitive count

  float cost = splitNodeCost + leftProb*nodeCost(leftTri,leftOther) + rightProb*nodeCost(rightTri,rightOther);
  //penalize duplication of primitives in both children (encourage disjoint splits)
  cost += duplicateCost*(leftTri + leftOther + rightTri + rightOther - totalTri - totalOther);
  //reward even splits for nodes with lots of triangles
  //see if max count reduced significantly below 3/4 of original
//  int reduction = totalTri+totalOther;
//  reduction = reduction>>1 + reduction>>2 - 64 - maxv(leftTri+leftOther,rightTri+rightOther);
  //clamp negative values to zero
//  reduction = reduction&~(reduction>>31);
//  cost -= reductionBenefit*reduction;
  return cost;
}
KDSplitCode KDBGreedyCostSplit::splitFixed(const AxisAlignedBox &bounds, const GeometrySet &geom)
{
  fixed.initialize(bounds);
  geom.estimateFixedSplits(fixed);
  fixed.integrateCounts();
  assert(fixed.getTotalTri()+fixed.getTotalOther()>=geom.size());
  float sizeX = bounds.max[0] - bounds.min[0];
  float sizeY = bounds.max[1] - bounds.min[1];
  float sizeZ = bounds.max[2] - bounds.min[2];
  assert(sizeX>0);
  assert(sizeY>0);
  assert(sizeZ>0);
  float area = (sizeX*sizeY + sizeX*sizeZ + sizeY*sizeZ);
  float denom = 1.0f/area;
  int totalTri = fixed.getTotalTri();
  int totalOther = fixed.getTotalOther();
  float bestCost = nodeCost(totalTri,totalOther);
  int bestIndex = 0;
  int bestType = KDST_LEAF;
  //test all possible recorded splits to find the best
  //starting with X splits
  float ax = (sizeY + sizeZ)*sizeX*denom;
  float bx = sizeY*sizeZ*denom;
  for(int i=1;i<KDFIXEDSPLIT_NUM;i++) {
    float relVal = relSplitTable[i];
    float cost = 
      splitCost( relVal*ax+bx,     fixed.getLeftTriCountX(i),fixed.getLeftOtherCountX(i),
                (1.0f-relVal)*ax+bx,fixed.getRightTriCountX(i),fixed.getRightOtherCountX(i),
                totalTri,totalOther);
    if (cost < bestCost) {
      bestCost = cost;
      bestType = KDST_SPLIT_X;
      bestIndex = i;
    }
  }
  //next do Y splits
  float ay = (sizeX + sizeZ)*sizeY*denom;
  float by = sizeX*sizeZ*denom;
  for(int i=1;i<KDFIXEDSPLIT_NUM;i++) {
    float relVal = relSplitTable[i];
    float cost = 
      splitCost( relVal*ay+by,     fixed.getLeftTriCountY(i),fixed.getLeftOtherCountY(i),
                (1.0f-relVal)*ay+by,fixed.getRightTriCountY(i),fixed.getRightOtherCountY(i),
                totalTri,totalOther);
    if (cost < bestCost) {
      bestCost = cost;
      bestType = KDST_SPLIT_Y;
      bestIndex = i;
    }
  }
  //next do Z splits
  float az = (sizeY + sizeX)*sizeZ*denom;
  float bz = sizeY*sizeX*denom;
  for(int i=1;i<KDFIXEDSPLIT_NUM;i++) {
    float relVal = relSplitTable[i];
    float cost = 
      splitCost( relVal*az+bz,     fixed.getLeftTriCountZ(i),fixed.getLeftOtherCountZ(i),
                (1.0f-relVal)*az+bz,fixed.getRightTriCountZ(i),fixed.getRightOtherCountZ(i),
                totalTri,totalOther);
    if (cost < bestCost) {
      bestCost = cost;
      bestType = KDST_SPLIT_Z;
      bestIndex = i;
    }
  }
  buildPredLeftTri = fixed.getLeftTriCount(bestType,bestIndex);
  buildPredRightTri = fixed.getRightTriCount(bestType,bestIndex);
  buildPredLeftOther = fixed.getLeftOtherCount(bestType,bestIndex);
  buildPredRightOther = fixed.getRightOtherCount(bestType,bestIndex);
  return fixed.getBestSplitCode(bestType,bestIndex);
}

KDSplitCode KDBGreedyCostSplit::splitAll(const AxisAlignedBox &bounds, const GeometrySet &geom)
{
  all.initialize(bounds);
  geom.estimateAllSplits(all);
  assert(all.getTotalTri()+all.getTotalOther()>=geom.size());
  float sizeX = bounds.max[0] - bounds.min[0];
  float sizeY = bounds.max[1] - bounds.min[1];
  float sizeZ = bounds.max[2] - bounds.min[2];
  assert(sizeX>0);
  assert(sizeY>0);
  assert(sizeZ>0);
  float area = (sizeX*sizeY + sizeX*sizeZ + sizeY*sizeZ);
  float denom = 1.0f/area;
  int totalTri = all.getTotalTri();
  int totalOther = all.getTotalOther();
  float bestCost = nodeCost(totalTri,totalOther);
  float bestSplitValue = 0;
  int bestType = KDST_LEAF;
  //test all possible recorded splits to find the best
  //starting with X splits
  float ax = (sizeY + sizeZ)*denom;
  float bx = sizeY*sizeZ*denom;
  int count = all.prepareXEvents();
  for(int i=0;i<count;i++) {
    float value = all.event[i].value;
    float cost = 
      splitCost((value-bounds.min[0])*ax+bx,all.event[i].leftTri,all.event[i].leftOther,
                (bounds.max[0]-value)*ax+bx,all.event[i].rightTri,all.event[i].rightOther,
                totalTri,totalOther);
    if (cost < bestCost) {
      buildPredLeftTri = all.event[i].leftTri;
      buildPredRightTri = all.event[i].rightTri;
      buildPredLeftOther = all.event[i].leftOther;
      buildPredRightOther = all.event[i].rightOther;
      bestCost = cost;
      bestType = KDST_SPLIT_X;
      bestSplitValue = value;
    }
  }
  //next do Y splits
  ax = (sizeX + sizeZ)*denom;
  bx = sizeX*sizeZ*denom;
  count = all.prepareYEvents();
  for(int i=0;i<count;i++) {
    float value = all.event[i].value;
    float cost = 
      splitCost((value-bounds.min[1])*ax+bx,all.event[i].leftTri,all.event[i].leftOther,
                (bounds.max[1]-value)*ax+bx,all.event[i].rightTri,all.event[i].rightOther,
                totalTri,totalOther);
    if (cost < bestCost) {
      buildPredLeftTri = all.event[i].leftTri;
      buildPredRightTri = all.event[i].rightTri;
      buildPredLeftOther = all.event[i].leftOther;
      buildPredRightOther = all.event[i].rightOther;
      bestCost = cost;
      bestType = KDST_SPLIT_Y;
      bestSplitValue = value;
    }
  }
  //next do Z splits
  ax = (sizeX + sizeY)*denom;
  bx = sizeX*sizeY*denom;
  count = all.prepareZEvents();
  for(int i=0;i<count;i++) {
    float value = all.event[i].value;
    float cost = 
      splitCost((value-bounds.min[2])*ax+bx,all.event[i].leftTri,all.event[i].leftOther,
                (bounds.max[2]-value)*ax+bx,all.event[i].rightTri,all.event[i].rightOther,
                totalTri,totalOther);
    if (cost < bestCost) {
      buildPredLeftTri = all.event[i].leftTri;
      buildPredRightTri = all.event[i].rightTri;
      buildPredLeftOther = all.event[i].leftOther;
      buildPredRightOther = all.event[i].rightOther;
      bestCost = cost;
      bestType = KDST_SPLIT_Z;
      bestSplitValue = value;
    }
  }
  if (bestType==KDST_LEAF) return KDSplitCode(KDST_LEAF);
  return KDSplitCode(bestType,bestSplitValue);
}

KDSplitCode KDBGreedyCostSplit::split(int depth, const AxisAlignedBox &bounds, const GeometrySet &geom)
{
  int size = geom.size();
  if (size==0) {
    return KDSplitCode(KDST_LEAF); //never split an empty node
  } else if (size <= KDALLSPLIT_MAX) {
    return splitAll(bounds,geom);
  } else {
    return splitFixed(bounds,geom);
  }
}


